  setwd("C:/Users/Jef/Desktop/Applied Statistics")
# setwd("D:/Mijn Documenten/applied statistics - project")
library(qcc)

#PROCESS CAPABILITY ANALYSIS:

paper <- read.table("paperfranklindata1.txt")

#informal check for outliers, Box and Whisker Plot
boxplot(paper, horizontal = TRUE, main="Box-and-Whisker plot of initial data")

#Kernel density plot
plot(density(paper$V1),main="Kernel density estimate of initial data",col="red",lwd=3)

#Normal probability plot
qqnorm(paper$V1,main="Normal probability plot of initial data", pch=19,cex=1,fg="red")
qqline(paper$V1,lwd=3,col="blue",lty="dashed")

#Confidence levels
samplenumber <- rep(1:35, each=1)
paperdata <- qcc.groups(paper$V1, samplenumber)
paperqcc <- qcc(paperdata,main="Process capability analysis of initial data", type="xbar",spec.limits=c(4.920, 4.980)) 
process.capability(paperqcc, spec.limits=c(4.920, 4.980))

#Shapiro-test
shapiro.test(paper$V1)

# proportion of non-conforming items:
# 1- P(LSL < X < USL) = \Phi(-3(2C_p - C_{pk})) + \Phi(-3C_{pk})
pnorm(-3*(2*1.195 - 1.02)) + pnorm(-3*1.02)

#CONTROL CHARTS
paper2 <- read.table("paperfranklindataallekollommen.txt")

#informal check for outliers, Box and Whisker Plot
boxplot(paper2[1:15,7], horizontal = TRUE, main="Box-and-Whisker plot of phase I data")


#Kernel density plot
plot(density(paper2[1:15,7]),main="Kernel density estimate of phase I data",col="red",lwd=3)

#Normal probability plot
qqnorm(paper2[1:15,7],main="Normal probability plot of phase I data", pch=19,cex=1,fg="red")
qqline(paper2[1:15,7],lwd=3,col="blue",lty="dashed")

#Confidence levels
paper2qcc <- qcc(paper2[1:15,2:6],main="Process capability analysis of phase1 data", type="xbar",spec.limits=c(4.920, 4.980)) 
process.capability(paper2qcc, spec.limits=c(4.920, 4.980))

shapiro.test(paper2[1:15,7])

#get center of first 15 samples
stats.R(paper2[1:15,2:6])
stats.xbar(paper2[1:15,2:6])

#get standard deviations of firs 15 samples
sd.xbar(paper2[1:15,2:6])
sd.R(paper2[1:15,2:6])

#calculating standard control limits x-bar chart
n	<-	5
mu	<-	4.958133
ER	<-	0.02566667
sigma	<-	0.01162417
D001	<-	0.367
D999	<-	5.484
d2	<-	2.326
LCLXs <-	mu-(3*sigma)/(sqrt(n))
UCLXs <- 	mu+(3*sigma)/(sqrt(n))
LCLRs <-	D001*ER/d2
UCLRs <- 	D999*ER/d2

#calculating the control limits they used
A2	<-	0.557
D3	<-	0.000
D4	<-	2.115
LCLXa <-	mu-A2*ER
UCLXa <-	mu+A2*ER
LCLRa <- 	D3*ER
UCLRa <- 	D4*ER

#getting the control charts with the limits used by the paper 
paperqccX1 <- qcc(paper2[1:15,2:6], type="xbar", center=mu, limits=c(LCLXa,UCLXa))
paperqccR1 <- qcc(paper2[1:15,2:6], type="R", center=ER, limits=c(LCLRa,UCLRa))

paperqccX2 <- qcc(paper2[16:32,2:6], type="xbar", center=mu, limits=c(LCLXa,UCLXa))
paperqccR2 <- qcc(paper2[16:32,2:6], type="R", center=ER, limits=c(LCLRa,UCLRa))

paperqccX3 <- qcc(paper2[33:47,2:6], type="xbar", center=mu, limits=c(LCLXa,UCLXa))
paperqccR3 <- qcc(paper2[33:47,2:6], type="R", center=ER, limits=c(LCLRa,UCLRa))


#CAPABILITY STUDY LAST 15 SAMPLES
paperqccnormaal <- qcc(paper2[33:47,2:6], type="xbar")
process.capability(paperqccnormaal, spec.limits=c(4.920, 4.980))